# Figure 1: generate the n-point example found in the paper
set.seed(102112516) 
N <- 50

X <- cbind(rnorm(N,0,1),rnorm(N,0,1)) # generate bivariate standard normal
X <- scaling.min.max(X)
ot <- otm2D(X,rep(1/N,N)) # obtain vector quantile 

# Calculate power diagram associated to "ot" to plot
pwd <- power_diagram(X[,1],X[,2],ot$weights)
centers <- t(sapply(1:N,function(i){colMeans(pwd$cells[[i]])})) # get cell centers

# Figure 1a: Vector Quantile
png('fig1a.png',width=8,height=8,units="in",res=300)
plot(c(0,1),c(0,1), xlab = expression("r"[1]),ylab = expression("r"[2]),type = 'n',
     cex.axis = 1.5, cex.lab = 1.5)
plot(pwd, add = TRUE)
for (i in 1:N){
  lines(rbind(X[i,],centers[i,]),lty = 2)
}
dev.off()


# Figure 1b: Lorenz Map Diagram
r <- c(0.6,0.4)
line1 <- rbind(c(r[1],0),c(r[1],r[2]))
line2 <- rbind(c(0,r[2]),c(r[1],r[2]))
A <- rbind(c(0.001,0.001),c(r[1],0.001),c(r[1],r[2]),c(0.001,r[2]),c(0.001,0.001))
box <-  st_polygon(list(A)) 

png('fig1b.png',width=8,height=8,units="in",res=300)
plot(c(0,1),c(0,1), xlab = expression("r"[1]),ylab = expression("r"[2]),type = 'n',
     cex.axis = 1.5, cex.lab = 1.5)
plot(box, col = rgb(1, 0, 0,0.25), add = TRUE)
for (i in 1:N){
 lines(rbind(pwd$cells[[i]],pwd$cells[[i]][1,]),col = 4,lwd = 1.5)
}
dev.off()